
#ifndef FASTSQRT_H_
#define FASTSQRT_H_

#include "libmath.h"

namespace p2 {

//#define _USE_FAST_MATH

#if defined(_USE_FAST_MATH)

/*
 * REF: http://en.wikipedia.org/wiki/Methods_of_computing_square_roots
 */
inline real32 fastSqrt(real32 z)
{
	union
	{
		int tmp;
		float f;
	} u;

	u.f = z;

	/*
	 * To justify the following code, prove that
	 *
	 * ((((val_int / 2^m) - b) / 2) + b) * 2^m = ((val_int - 2^m) / 2) + ((b + 1) / 2) * 2^m)
	 *
	 * where
	 *
	 * val_int = u.tmp
	 * b = exponent bias
	 * m = number of mantissa bits
	 *
	 * .
	 */

	u.tmp -= 1 << 23; /* Subtract 2^m. */
	u.tmp >>= 1; /* Divide by 2. */
	u.tmp += 1 << 29; /* Add ((b + 1) / 2) * 2^m. */

	return u.f;
}
/*
 * REF: http://en.wikipedia.org/wiki/Fast_inverse_square_root
 */
inline real32 fastInvSqrt(real32 x)
{
	real32 xhalf = 0.5f * x;
	uint32 i = *((uint32*) &x);
	i = 0x5f375a86 - (i >> 1);
	x = *((real32*) &i);
	return x * (1.5f - xhalf * x * x); // repeat this statement for a better approximation
}
inline real64 fastInvSqrt(real64 x)
{
	real64 xhalf = 0.5 * x;
	uint64 i = *((uint64*) &x);
	i = 0x5fe6eb50c7aa19f9 - (i >> 1);
	x = *((real64*) &i);
	x = x * (1.5 - xhalf * x * x); // repeat this statement for a better approximation
	return x;
}
#	ifdef _MSC_VER
//----------------------------------------------------------------------------
// REF: http://www.codeproject.com/Articles/69941/Best-Square-Root-Method-Algorithm-Function-Precisi
inline real64 fastSqrt(real64 x)
{
	__asm
	{
		fld x
		fsqrt
	}
}
inline real32 fastSqrt(real32 x)
{
	__asm
	{
		fld x
		fsqrt
	}
}
inline real64 __declspec (naked) __fastcall fastSqrt1(real64 x)
{
	_asm fld qword ptr [esp+4]
	_asm fsqrt
	_asm ret 8
}
#	endif
#endif

}

#endif	//FASTSQRT_H_
